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This  paper  deals  with  the  problem  of  optimization  of  the  final  turn-into-the-wind  maneuver  of  an 
aerial  delivery  system  with  account  of  the  best  known  winds.  The  wind  model  required  for  the 
optimization  algorithm  to  work  may  utilize  onboard  wind  estimates  only,  incorporate  the  ground 
winds  provided  a  priori  or  on-line  by  the  target  ground  station,  or  be  based  on  the  winds  measured  and 
uplinked  by  the  preceding  system.  The  previous  work  by  the  authors  took  care  of  the  major  touchdown 
error  contributor,  downwind  variation  of  the  winds.  The  effect  of  these  variations  was  mitigated  by 
constantly  recomputing  an  optimal  reference  trajectory  to  complete  a  final  turn  in  a  given  time.  This 
paper  presents  some  modifications  of  the  original  optimization  routine  to  accommodate  some  specific 
applications  including  intentional  landing  with  a  substantial  crosswind  component  and  operating  in  the 
mountainous  areas  with  significant  variations  in  the  vertical  component  of  the  wind  (updrafts  and 
downdrafts).  Specifically,  the  paper  presents  derivation  of  equations  to  account  for  one-,  two-  and 
three-dimensional  structure  of  the  wind.  In  addition,  adjustments  to  the  optimal  control  problem  using 
the  direct-method-based  approach  developed  earlier  for  a  simple  one-dimensional  wind  model  are 
developed. 


Abbreviations 


ADS 

= 

Aerial  Delivery  System 

AGL 

above  the  ground  level 

BCs 

= 

boundary  conditions 

CEP 

= 

Circular  Error  Probable 

FAC 

= 

final  approach  capture  (point) 

GPS 

= 

Global  Positioning  System 

TI 

= 

turn  initiation  (point) 

TPBVP 

= 

two-point  boundary- value  problem 

I.  Background 

IN  an  attempt  to  mitigate  the  effect  of  unknown  variable  winds,  Siegers  and  Yakimenko  formulated  the  following 
two-point  boundary- value  problem  (TPBVP)  (Fig.  1). 1,2  Using  a  right-handed  coordinate  system  {W}  aligned  with 
the  prevailing  ground  wind  (defining  a  downrange  axis)  we  need  to  bring  a  non-powered  aerial  delivery  system 
(ADS)  from  some  initial  point,  with  the  state  vector  defined  at  t  -  0  as 

xo  =[Wo>n]  (1) 

(x  -  downrange, y  -  crossrange  and  y/-  heading  in  {W})  to  another  point 

(2) 

at  t  -  tf  .  In  Eq.(2),  V*h  is  the  estimate  of  a  horizontal  component  of  a  steady-state  airspeed,  W  =  const  is  the  only 
component  of  the  wind  vector 

w  =  [r,0,0f  (3) 
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and  I*™  is  the  desired  final  approach  time. 

Figure  1  shows  the  portion  of  a  guided  descent  to  be  optimized  (appearing  in  between  two  vertical  lines).  It 
occurs  between  the  turn  initiation  (TI)  point  at  some  altitude  h0  (defined  by  the  estimate  of  W,  and  the  components 
of  the  ADS  velocity  vector  as  explained  in  Ref.l)  and  final  approach  capture  (FAC)  point. 


Figure  1.  Final  turn-into-the-wind  maneuver  and  landing  approach. 


The  altitude  at  t  =  tf  was  defined  based  on  the  constant  descent  rate  assumption 

\=vX;=K-tfv;  (4) 

(here  V*  is  the  estimate  of  a  steady-state  descent  rate). 

Hence,  we  need  to  find  the  trajectory  that  satisfies  the  boundary  conditions  (BCs)  (1)  and  (2)  along  with  the 
constraint  imposed  on  the  control  (turn  rate),  \y/\  <  ^max  ,  and  allows  completing  the  maneuver  in  exactly 

h 

T  =t  =-^-Tdes  (5^ 

1tum  V  1app  yJ) 

V  v 

The  assumption  of  a  constant  descent  rate  allows  eliminating  the  differential  equation  for  an  altitude  and 
reducing  ADS’  kinematics  down  to 

r  V  r*r»Q  nr  W 

(6) 

From  these  two  equations  it  follows  that  if  the  final-turn  trajectory  is  given  (defined  analytically  by  x(t )  and 
y(t ) ),  then  the  yaw  angle  along  this  trajectory  is  related  to  the  change  of  the  inertial  coordinates  as 


X 

Vl  cos  y/  +  W 

_y_ 

V*  sin  i// 

y/  -  tan 


y 

x-W 


(7) 


Differentiating  Eq.(7)  provides  with  the  yaw  rate  control  required  to  follow  the  reference  fmal-tum  trajectory  in  a 
presence  of  a  constant  downwind  W 


y(x-W)-xy 
(x-W)2 +y2 


(8) 


Once  the  TPBVP  is  solved,  Eq.(8)  provides  the  time  history  of  the  optimal  control  y/opt  (t)  .  It  is  this  control 

profile  that  is  tracked  by  the  ADS’  control  unit  as  explained  in  Ref.l. 

The  developed  guidance  and  control  algorithm  was  implemented  on  the  Snowflake  ADS.3  In  the  period  between 
May  of  2008  and  May  of  2011  this  system  has  been  dropped  from  different  deployment  platforms  from  altitudes 
2,000-14,000  ft  above  the  ground  level  (AGL)  over  150  times.2  During  the  first  set  of  three  drops  in  May  of  2008 
the  Snowflake  ADS  exhibited  the  circular  error  probable  (CEP)  of  55m  with  the  standard  deviation  of  9m.4  These 
parameters  were  gradually  reduced  to  the  CEP  of  1  lm  with  the  standard  deviation  of  6m,  exhibited  in  the  set  of  four 
drops  in  August  of  20 10. 4  This  outstanding  performance  of  the  smallest  autonomously  guided  ADS,  featuring  the 
cheapest  and  therefore  the  worst  canopy  and  being  most  susceptible  to  the  winds  was  achieved  by  implementing  the 
latest  technologies  in  control  theory1,6  and  also  by  utilizing  the  best  available  options  of  accounting  for  the  unknown 
surface-layer  winds.3,7,8 
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This  paper  presents  modifications  of  the  original  optimal  control  (8)  based  on  alternative  wind  models,  which 
can  be  quite  different  from  the  one  presented  in  Eq.(3).  To  this  end,  the  following  section  presents  some  flight  test 
data  recorded  by  an  onboard  sensors  suite  during  a  couple  of  arbitrary  chosen  drops  showing  that  the  actual  values 
of  V*  and  W  can  vary  quite  drastically.  Based  on  observations  of  these  data,  Sections  III  and  IV  discuss  different 

application-specific  modifications  of  the  final-turn  optimization  routine  of  Ref.l  and  briefly  described  in  the 
beginning  of  this  section.  More  precisely,  Section  III  accommodates  linear  and  logarithmic  wind  profiles  in  the 
downwind  direction  (ID  winds  optimization),  and  Section  IV  considers  both  downwind  and  crosswind  components 
(2D  optimization).  Section  V  simplifies  computational  algorithms  of  Section  IV  by  utilizing  the  precomputed 
ballistic  winds.  Section  VI  addresses  variations  in  the  initial  conditions  for  the  final-turn  maneuver  caused  by 
variable  2D  winds.  Section  VII  presents  a  simple  way  of  mitigation  the  effects  of  wind  updrafts  and  downdrafts  (3D 
optimization). 


II.  Surface-Layer  Winds 

Figure  2  presents  two  samples  of  data  measured/estimated  and  recorded  onboard  the  Snowflake  ADS.  The  plots 
in  the  first  column  belong  to  one  drop  terminated  with  the  miss  distance  of  10m,  and  in  the  second  column  -  to  a 
second  drop  that  resulted  in  a  miss  distance  of  6m.  The  first  row  of  the  plots  presents  the  altitude  versus  time  profile. 
The  second  row  shows  estimated  downwind  component  of  the  wind  (which  changes  all  the  time  based  on  the  latest 
observations  of  the  ground  track  speed  measured  by  the  onboard  Global  Positioning  System  (GPS)  receiver).  The 
third  row  shows  vertical  speed  of  the  ADS  as  measured  by  GPS  and  smoothed  (filtered)  by  the  onboard  inertial 
navigation  unit.  The  last  row  of  the  plots  presents  different  stages  of  flight  when  these  data  were  collected. 


Figure  2.  Flight  parameters  recorded  during  the  two  drops  in  August  of  2010. 

As  seen  in  practice,  neither  Vv*  nor  W  used  in  equation  of  Section  1  are  constant.  While  the  first  set  of  data 
exhibits  a  kind  of  gradual  decrease  of  the  downwind  component  W  with  time  (altitude),  the  second  set  of  data 
features  more  or  less  constant  winds  up  to  about  200m  altitude  AGL  with  a  sudden  halved  decrease  below  this 
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altitude.  For  convenience  of  analysis,  Fig.3  presents  the  values  of  V*  and  W  for  the  same  sets  of  data  as  versus 

altitude,  rather  than  time.  Looking  at  both  sets  of  data  together  one  can  notice  that  the  descent  rate  (negative  of  the 
vertical  speed)  varies  from  Om/s  (as  a  result  of  some  updraft  motion)  all  way  up  to  6m/s,  apparently  affected  by  a 
close  to  the  ground  downdraft.  It  is  indicative  that  both  sets  of  data  belong  to  the  same  ADS  dropped  at  the  same 
location  less  than  an  hour  apart.  No  wonder  that  parafoils  systems,  being  influenced  by  varying  winds  exhibit 
inconstant  performance. 


Figure  3.  Altitude  dependences  of  the  ADS  estimates. 


According  to  the  guidance  strategy  described  in  Ref.  1  the  winds  aloft  only  affect  the  location  of  the  computed  TI 
point  along  the  downwind  leg.  This  leaves  all  unmodeled  dynamics  to  be  handled  at  the  following  stages  of  flight, 
i.e.  final  (base)  turn  and  final  landing  approach.  Data  presented  in  Figs.2  and  3  is  zoomed  to  the  final  stages  in 
Figs.4  and  5  to  show  parameter  variations  at  the  surface  layer,  staring  at  about  300m  to  include  the  downwind  leg 
where  the  decision  to  turn  is  made. 

Obviously,  in  the  general  case  the  unaccounted  winds  may  have  components  in  all  three  directions 

W  dls‘(h)  =  [wx,wy,wzJ  (9) 

(here  wx  denotes  a  downwind  component,  not  accounted  for  by  Eq.(3),  while  wz  is  considered  positive  for 

downdrafts  to  be  consistent  with  the  descent  rate  sign  convention).  With  disturbances  (9)  the  kinematic  equations  (6) 
become  three-dimensional 


X 

V*h  cos  y/  +  W  + 

y 

= 

Vlsmy/+wy 

h 

-K-k 

In  the  original  algorithm  of  Ref.  1  however  the  unaccounted  winds  (9)  were  treated  as  disturbances,  so  that  it  was 
up  to  the  control  system  to  mitigate  their  effect  while  still  using  Eq.(6)  to  compute  the  reference  control  (7).  Surely, 
these  disturbances  were  then  the  primary  reason  for  the  ADS  not  tracking  the  calculated  optimal-turn  trajectory 
precisely.  As  shown  in  numerous  simulations  and  in  practice  it  is  these  winds  that  can  cause  the  ADS  to  land  short 
of  the  target  (in  the  case  of  the  higher  than  expected  head  winds)  or  overshoot  it  (tail  winds).  That  is  why  the  optimal 
trajectory  needs  to  be  constantly  updated  during  the  final  turn,  each  time  starting  from  the  current  (off  the  original 
trajectory)  initial  conditions  (IC)  and  still  forcing  the  ADS  to  be  at  point  (3)  within  an  updated  Tturn  (5). 

The  goal  of  the  following  sections  is  to  account  for  wind  disturbances  (9)  at  the  stage  of  generating  the  reference 
control,  i.e.  trying  to  use  Eq.(10)  instead  of  Eq.(6).  Obviously,  it  can  be  done  only  if  the  wind  disturbance 
components  (9)  can  be  modeled  (using  more  information  about  the  winds  known  a  priori).  To  this  end,  Section  III 
starts  with  more  accurate  modeling  of  downwind  component  of  the  surface  winds,  followed  by  Section  IV 
introducing  a  crosswind  component  and  ending  with  Section  V  discussing  the  vertical  wind  component. 


III.  Optimization  Based  on  the  Linear  and  Logarithmic  Surface-Layer  Wind  Models 


Assume  that  instead  of  a  constant  v-component  of  the  prevailing  wind  W  versus  altitude  h  (Eq.(2))  we  have  a 
linear  profile 


W(h)  =  WG+Kh 


(11) 
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where  WG  is  a  known  ground  wind  and  coefficient  K  -  (W0  -WG)h0l  is  defined  by  the  ground  wind  WG  and  wind 
W0  measured  at  an  altitude  h0  (corresponding  to  the  point  where  the  final  turn  begins).  In  terms  of  Eq.(9)  it  means 
that  we  are  trying  to  model  the  downwind  disturbance  as  wx  ( h )  =  WG  -W0  +  Kh  .  Such  a  profile  might  be  based  on 

the  known  ground  winds  (available  from  the  nearby  airport,  measured  by  the  target  ground  station,3,7  etc.)  uplinked 
to  the  descending  ADS. 


b) 

Figure  4.  Flight  parameters  exhibited  at  the  300m  surface-layer  (zoomed-in  versions  of  Fig.2). 
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Figure  5.  Altitude  dependences  of  the  ADS’  estimates  between  the  surface  and  300m  AGL  altitude. 


Then,  the  original  TPBVP  of  Section  I  should  be  reformulated  for  a  slightly  different  system  of  kinematic 
equations 
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X  =  ^  cos  Y  +  W(h) 

_y\  _  Vl  sin  y/ 

and  different  BCs.  To  be  more  specific,  starting  from  some  (different)  initial  point,  defined  by  a  different  expression 
for  a  distance  past  the  target  on  the  downwind  leg  to  initiate  the  final  turn  maneuver  (defined  in  Ref.  1  as  Dswitch ), 
we  will  need  to  bring  a  parafoil  to  the  point 

x/  =[(K-w)t%;  -±kvv\t£;)2,o,-x]t  (i3) 

(computation  of  Dswitch  will  be  addressed  in  Section  VI). 

To  compute  the  offset  in  Eq.(13)  we  used  the  fact  that  the  final  landing  approach  starts  at  the  altitude  V*T^ ,  so 
that  using  an  obvious  relation 

t  dh 

dt- - -  (14) 


we  may  write 


x r  =  K -W  dt=  K  -W 


vapp  „ 

J  (vl-W)^  =  (vl-WG)T^;-\KV:{Tda;sp)2  (15) 


In  this  case,  inverting  equations  (12)  yields 


^=tan-' — ^ — 
x-W(h ) 


Compared  to  Eq.(7),  Eq.(8)  features  and  altitude-dependent  wind  profde  W(h) ,  so  that  its  differentiation  with 
account  of  Eq.(l  1)  results  in  a  slightly  different  equation  for  the  turn  rate 

y(x  -W)-(x-  W)y  y(x  -W)-(x-  KV*)y 


(. x-W)2+y 2 


(x-W)2  +y2 


(cf.  with  Eq.(8)). 

The  only  modifications  the  numerical  algorithm  described  in  Ref.  1  requires  in  this  case  is  that  it  should  involve 
the  new  BCs 


K  cos  t^o  +  ^o 

X 

-0 

V*hsm¥  o 

Wr=0 

*  =  (K-wG)T%;-tKv;(T%;)2 1  Tx  =  -i v;+wg+kvx 

-y\r=T,  [  0  J’  \y\r=Tf  [  0 

as  well  as  computation  of  an  altitude 

hj  =hj_1-V*Atj_l,  j  =  2,...,N  =h0) 

and  the  corresponding  wind  magnitude  at  each  computational  node 

Wj  =WG+  Khj 

The  latter  two  values  then  are  to  be  used  to  compute  time  intervals  between  two  computational  nodes 


AtJ~ 1  \  ts*2 


(xj  —  Xy_!  )2  +  (y j  —  y j-\  )2 


and  heading 


Vh z  +  0 25(Wj  +  Wj_x )  +Vh  (1 Wj  +  Wj_x ) cos ¥]_x 

...  _  hyj 


y/ .  =  tan 


Vrwi 


In  the  case  when  the  ground  winds  are  not  available  a  general  logarithmic  wind  profile  may  be  used  in  lieu  of  the 
liner  profile  of  Eq.(l  l)8 

W(h)  =  p  +  a\n{h)  (24) 

In  this  case  some  of  the  above  equations  will  be  replaced  with  the  new  ones 
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K  TZ 


Xf=  J  {<-W)y  =  {vl+a-p) T%;  - aT^  ln(Kv*7^ ) 


X 

_y_ 

T=  0 

-Wo  vh  sin  ^0-^o  K 
Vo K  cos^o 


\K in//p 


X 

_y_ 

T=Tf  _ 

-v; +p+a\n{VX;sp) 
0 


and 


Wj  =  ,6  + a  In/) 


(25) 

(26) 

(27) 

(28) 


The  remaining  equations  will  still  be  the  same. 


IV.  Accommodating  Cross-Wind  Data 

The  optimization  routine  of  the  Snowflake  guidance  algorithm  can  also  accommodate  crosswinds  if  they  are 
known  in  one  form  or  another.  This  capability  may  be  useful  in  organizing  a  swarm  attack  or  landing  onto  a  ship’s 
deck,  not  necessarily  aligned  with  the  wind.9  Consider 

wx(h)  =  fx(h ),  wy(h)  =  fy(h)  (29) 

to  be  x-  (downwind)  and  y-  (crosswind)  components  of  a  horizontal  wind  profile  approximated  with  some  analytical 
dependence  (e.g.  of  the  form  of  Eq.(l  1)  or  (24)  or  cubic  spline).  Then,  we  can  write 

(30) 

Note  that  instead  of  Eq.(6)  we  are  now  using  the  first  two  equations  of  (10),  emphasizing  that  the  entire  trajectory  is 
intentionally  aligned  not  with  the  major  wind  component,  so  that  w  (h)  can  actually  be  even  larger  than  wx{h)  (i.e. 

we  let  W  =  0 ). 

Accounting  for  the  new  kinematics  described  by  Eq.(30)  the  final  point  can  be  defined  as 

(31) 


X 

vl  cos  y/  +  wx  (h) 

_y_ 

sin  i//  +  wy(h)_ 

K  TZ 


where  the  offsets  in  x-  and  y-  direction  will  be  computed  as 

t r*  rjides 

V  aPP  TT 

T/*rpdeS  f  /7\  f  /  7  \  dfl 

xf=vhTaPP-  J  wx(h)~pr,  yf=-  J  wv(7i)— 
0  Vv  0  Vv 

The  heading  angle  equation  will  be 


y/  =  tan 


■  1  y-Wy(h) 
x-wjh) 


while  its  derivative  will  be  presented  by 


¥  =  - 


(y  -  w !y  (h)Vv  )(i  -  wx  (h))  -(x-w'x  (h)Vv  )( y  -  wy  (/»)) 
(x-wx(h))2+(y-wJh))2 


(where  Wx  =  dwx  /  dh  and  wy  =  dwy  /  dh).  The  total  speed  will  now  be  expressed  as 

I  vg  I  =  V*2  +  y1  =  2  +  ( Vc  W +^wf  +  zK  ( Vc  W cos  v+Wy  W sin  v) 


The  numerical  procedure  will  proceed  with  the  boundary  conditions 


X 

xo 

X 

_y_ 

T= 0 

Xo_ 

? 

_y_ 

T=0 

Vh  cos^o  +wx(h0) 

K  sin  Vo  +  H’v)/;o  ) 


r= 0 


~Vo vh  sin^o  +K(h)Vy 
Vo  K  cos^q  +  w'y(h)V* 


X 

V 

X 

-K +wx(h ) 

X 

"0" 

_y_ 

T=Tf 

y/. 

? 

_y_ 

T=Tf 

Wy(tl) 

? 

_y_ 

T=Tf 

0 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 
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(38) 


and  will  involve  computing  wind  components  at  each  step 

Wxj  =  fx  (hj  )  and  =  fy  (hj  ) 

to  be  used  in  the  numerical  equations  similar  to  those  given  in  Section  IV. 


V.  Using  Ballistic  Winds 

In  the  previous  section  we  were  relying  on  some  analytical  wind  profiles  wx  =  fx(h)  and  w  =  fy(h)  .  In 
practice  however  the  components  of  the  horizontal  wind  can  be  available  as  the  look-up  tables,  containing  triples  hj  , 
wxj  ,  and  wyj  (or  wind  magnitude  and  direction,  instead  of  the  last  two  parameters).  In  this  case,  to  avoid  computing 
derivatives  wx  and  wf  in  Eqs.  (34)  and  (36)  we  could  use  the  so-called  ballistic  winds.  By  definition,  if  at  some 
altitude  //we  have  a  ballistic  wind  of  magnitude  WB  and  direction  x¥w  ,  then  the  effect  of  variable  winds  wx(h)  and 
w  (h)  for  some  system  with  the  descent  rate  Vv  on  its  way  from  altitude  H  down  to  the  surface  is  reduced  to  simple 
formulas 


or  in  other  words 


m  =  yWB  COS4V,  y(h)  =  yWBsm'Vw 


J wx (h)dh  =  ^WB  cos  ^ ,  ]^Wy(h)dh  =  ^ WB  sin  ¥ 


0K 


0K 


V 


(39) 


(40) 


(Note,  usually  wx  and  wy  are  measured  in  the  local  North-East-Down  coordinate  system,  so  ¥w  is  a  direction  with 
respect  to  the  true  North.  However  in  our  case  wx  and  w  are  expected  to  be  provided  by  the  first  descending 
system  in  {W},  so  that  Ww  will  also  be  calculated  in  {W}). 

Substituting  definite  integrals  in  Eq.(40)  with  the  finite  sum  of  trapezoids  based  on  the  discrete  values  of  hk  ,  wk 
and  y/wk,  k  =  1,...,M  ,  we  get 


£(A*-V.) 


Wj*  +  w  cos\p  V  (h  -h  ) Wyk  -  -  h  w  sin^p 

0  nMVV  M  X  W;M  ■>  Z^Vlk  nk-l  )  0  nMVV  M  'T  u 


(41) 


,LMrr  M  0111  A  W;M 

k= 2  ^  k= 2  ^ 

The  index  starts  from  2  because  by  definition  the  winds  measurements  at  the  lowest  altitude  can  be  considered 
ballistic  winds  at  this  altitude. 

From  Eqs.(41)  it  further  follows  that 


tan'E 


M 

I 

k= 2 
W;M  M 


YJ{K-K-l)[Wyk+Wy,k-l) 


Yk(hk-hk-l){Wxk+Wx,k-l) 


(42) 


f  M 


=- 


2K 


\k= 2 


f  M 


\k= 2 


For  the  specific  case  when  hk  -  hk_x  =  Ah  =  const ,  k  -  2 ,...,M  ,  Eqs. (42)  can  be  further  reduced  to 


tan'E 


M 

ZE,*+wM-i) 


W;M  M 


w  = 

VV  M 


Ah 


£(wxt+wxfi-l) 


2  h 


M  M 

+  Wx,k- 1  )  +  Z  (Wyk  +  Wy,k- 1  ) 


lM  VV*=2 


(43) 


If  the  ballistic  winds  are  known  a  priori ,  meaning  that  wx  =  fx(h)  and  wy  =  f  (h )  were  provided  by  the  first 

descended  system  as  a  look-up  table,  then  the  original  guidance  algorithm  of  Ref.  1  assuming  constant  winds,  can  be 
used  with  no  change. 
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VI.  Computation  of  the  Initial  Conditions  for  the  Final  into-the  Wind  Turn 

In  addition  to  the  modifications  of  the  original  algorithm  of  Ref.  1  discussed  in  the  previous  section,  different 
wind  models  require  changing  the  initial  conditions  to  initiate  the  final  turn  as  well.  Following  the  original  guidance 
algorithm  of  Ref.l,  the  x-  (downwind  direction)  budget  equation  for  two  phases,  base  turn  (t  e[t0;tx])  and  final 

approach  (t  ^\tyt2])  can  be  represented  as 


ll  l2  l2 

Dswitch  =  -J  W(h)dt  +  \(V*-W(h))dt  =  V*hTapP~\w(h)dt  =  V*hTapp -(Ttum  +  Tapp)W 


-.|0 

K 


(44) 


k)  k  k) 

Here  Dswitch  is  the  distance  passed  the  target’s  traverse  when  the  base  turn  should  be  initiated  (in  this  case  upon 

completion  of  the  aforementioned  two  phases  the  ADS  will  be  right  on/above  the  target),  and  W 1°  denotes  a 

% 

downwind  component  of  the  wind  averaged  within  the  altitude  range  h  e  [0;/z0] . 

The  altitude  budget  equation  for  these  two  plus  the  portion  of  the  downwind  leg  phase  starting  at  some  altitude  h 
at  a  distance  x  from  the  target’s  traverse  is 


f 


h  =  V, 


-x  +  A 


switch 


+  y  T  +V  T 

r  v  turn  v  app 


(45) 


h  J 


Ideally,  when  x  =  Dswitch  the  altitude  h  should  be  equal  to  V*  ( Tturn  +  Tapp ) .  In  practice  however  it  may  not  be  the 
case.  In  order  to  eliminate  the  altitude  error,  we  may  want  to  adjust  the  actual  final  approach  time  Tapp  . 

Resolving  Eqs.  (45)  and  (45)  with  respect  to  Dswitch  and  Tapp  yields 


D, 


switch 


)  -xir: + w 


S«|0 

K 


+  TturnVh\Vh-W  | 


^|0  -.0 


K 


^|0  -.0 


V„nrh-2W\K  + 


w[ 


2Vh  -2WL  +W | 


•  o 
K 


,  .  ^  —  o  —  o 

h(Vh-W  +W 


—x  +  Tt„ 


*  — 10  —  0 


app 


0  -  0  \ 

+  W 

2V*h-2  W°  +W{ 

K  h) 

h  K  1 

(46) 


(47) 


(here  we  used  an  obvious  relation  Wj^°  =  W\h  -  ). 

In  the  case  of  W  -  W 1°  =  w\°  =  const(h ) ,  Eqs.  (46)  and  (47)  are  simplified  to  those  of  the  original  guidance 

\h  \h0 

algorithm  of  Ref.l.  In  all  other  cases,  Eqs.  (46)  and  (47)  have  to  be  used.  For  the  linear  and  logarithmic  surface-layer 
wind  models  (12)  and  (24)  of  Section  III  the  averaged  winds  from  some  current  altitude  h  down  to  the  ground  can  be 
computed  as 


w\°  = 


1 

w\°h  =-j(WG+  Kh )  dh  =  WG  +  j  Kh 

0 

h 

-  +  a\n(h))dh  =  /3  +  a\n{h)-a 


(48) 


(49) 


respectively.  Substituting  JV\°  and  w\°  ,  computed  using  Eq.(48)  or  Eq.(49),  into  Eqs. (46), (47)  results  in  in  a  wind- 

\h  \h0 

model-specific  values  for  Dswitch  and  Tapp  (they  are  quite  bulky  and  are  not  given  here). 

Alternatively,  the  values  for  W 1°  and  W 1°  can  be  substituted  with  the  corresponding  downrange  component  of 

\h  % 

the  ballistic  winds  (computed  in  accordance  with  Eq.(42)  for  h  and  h0  ,  respectively) 
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wl=WB(h)cos(Vw{hj),  W^=WB(h0)Cos('Vw(h0))  (50) 

In  the  case  of  the  two-dimensional  wind  model  (including  the  crosswind  component),  Eqs.(46)-(50)  will  still  be 
valid,  since  computation  of  Dswitch  relies  on  the  downrange  winds  only. 


VII.  Accounting  for  Vertical  Wind  Disturbances 

Suppose  we  managed  to  have  two  ADSs,  so  that  while  descending  the  first  one  produces  and  passes  the  winds 
estimates  to  the  second  one.  Such  estimates  will  be  represented  by  the  GPS  time- stamped  quadruples:  hk ,  wxk ,  w  k  , 

and  wzk,  k-  2,3 ,...,M  (k  =  1  corresponds  to  h  =  0).  As  shown  in  the  previous  sections,  even  if  triplets  hk ,  wxk , 
wyk  are  available,  accounting  for  these  data  while  generating  a  reference  trajectory  may  still  pose  a  computational 
problem.  Accounting  for  the  vertical  component  of  the  wind  wz(h) ,  i.e.  updrafts  and  downdrafts,  is  even  more 
complicated.  In  a  non-stable  atmosphere  this  component  of  the  wind  may  cause  the  same  type  of  problem  as 
unaccounted  for  horizontal  winds  simply  because  it  changes  the  descent  time  forcing  ADS  to  land  sooner  (shorter  of 
the  target)  and  later  (resulting  in  the  overshoot). 

For  example,  consider  a  sudden  updraft  on  the  downwind  leg  (Figs.4a  and  5a)  or  downdraft  while  ADS  is  at  the 
final  turn  (Figs.4b  and  5b).  Obviously,  such  events  may  cause  a  serious  problem.  At  the  downwind  leg  a  vertical 
motion  of  the  air  mass  cause  a  violation  of  the  altitude  budget  equation  (45).  The  final  turn  maneuver  is  therefore  the 
last  chance  to  mitigate  this  violation.  However,  updrafts  and  downdrafts  at  this  phase  of  descent  mess  the  optimal 
solution  obtained  assuming  a  certain  time  of  maneuver  Tturn  (Eq.(5)).  The  capability  to  recompute  this  maneuver 
while  turning  even  at  each  control  cycle  (if  needed)  was  a  real  breakthrough  of  the  original  guidance  algorithm.1  Yet, 
as  opposed  to  updrafts,  downdrafts  may  be  still  of  the  major  concern  because  they  may  decrease  the  time  of  the 
fmal-tum  maneuver  at  once,  leaving  no  time  to  recover.  Hence,  accounting  to  the  vertical  winds  may  be  quite 
beneficial. 

Suppose  that  the  vertical  component  of  the  wind  is  known.  Again,  it  most  likely  comes  as  a  look-up  table,  hk  vs. 
wzk ,  but  theoretically  we  could  use  a  low-order  polynomial  regression  to  approximate  it  with  some  analytical 
dependence  wz{h)  .  In  this  case  this  dependence  may  be  used  to  modify  the  vertical  motion  equation  (14)  to 


dt  =  - 


dh 


(51) 


Vv  +wz(h) 

This  equation  is  then  to  be  used  in  Eqs.  (15),  (25),  (32),  (40),  (45),  (48),  and  (49).  Obviously,  depending  on  the 
specific  analytical  representation  wz  ( h )  the  resulting  equations  may  be  very  bulky,  so  the  alternative  approach  may 
be  based  on  the  analogous  of  the  ballistic  winds  concept  introduced  in  Section  V,  which  does  not  require  analytical 
regression  but  can  rather  utilize  ( hk ,  wzk )  pairs  explicitly. 

Following  Eqs.  (40)  and  (41),  let  us  introduce 

1  r  1  M 

=  -  j  wz(h)dh  ~  )(HA  +  w-k-\ )  (52) 

0  k=2 

which  denotes  a  downdraft  component  of  the  wind  averaged  within  the  altitude  range  h  e  [0 ;h\ .  Using  this  average 
downdraft  we  can  rewrite  Eq.(51)  as 

(53) 


This  equation  (implicit  in  h )  allows  you  to  estimate  time  T  needed  to  descent  from  the  altitude  h.  Now  let  us  use 
Eq.(53)  to  correct  Eq.(5).  To  this  end  let  us  use  Eq.(53)  to  replace  Eq.(4)  with 


/  T7-*  -  0 

v v  +wzL 


\fdes  _  h 
l)1app  n0 


Noting  that  wz  hp' 


_  o 

:  Wz  , 
hTi 


_  0 
-Wz  . 

z  hF, 


we  arrive  to 


ZK, 


h0-( V*  +  w. '° 
T  = - - 

turn  (  * 

(r- 


—  0 

+  W^r, 


_  0 
Wz 

z  hFf 


(54) 


(55) 
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Using  this  corrected  value  of  the  final  turn  maneuver,  allows  onboard  guidance  unit  to  produce  a  more  balanced 
control  input  \jr  (/) . 

In  addition  to  correcting  the  time  of  the  maneuver  the  optimization  algorithm  may  include  an  additional  varied 
parameter.  To  be  more  specific,  we  can  replace  the  final  condition  (2)  with  the  new  one 

X/  =  |>;  -  W)T%  cos  ¥f ,  -{v;  -  W)T£  sin  ¥f ,  ¥f  J  (56) 

with  the  value  of  y/f  being  variable.  This  allows  cutting  the  final  turn  maneuver  and  gliding  directly  to  the  target 
with  y/f  ^  -7i .  Similar  changes  can  be  made  for  any  wind  model  considered  in  the  previous  sections.  Specifically, 
the  final  conditions  of  Eq.(19)  can  be  substituted  with 


X 

-\kv:{t^sp)2  y0Wf 

X 

A  cos ¥f+wG+Kvv*Tad;;~ 

_y_ 

T=Tf 

-((K-wg)P;;  Ak<(t^)2)  sm¥f 

? 

_y_ 

T=Tf 

f  sin  Y  f 

Eq.(27)  with 


X 

_y_ 

T=Tf 

and  Eq.(37)  with 


+  a-p)T%sp  - ccT^pp  \n(V* )) cos ¥ f 
,+a-p)Tl -aT%;  ln^Xpjsin^ 


X 

2  _j_  yl  CQS^ 

_y_ 

T=Tf 

-f2f+y2fsm¥f_ 

X 

K  cos  Vf  +  P  +  a  HKTw  ) 

? 

_y_ 

T=Tf 

L  f  sin  ¥ f 

J 

1  _ 

Vl  cos  y/f  +  wx  (A)l 

V  L 


V*hsmy/f+wy(K) 


,  (58) 


(59) 


VIII.  Conclusion 

The  direct-method-based  approach  to  optimize  the  final-turn  maneuver  and  mitigate  all  the  errors  of  the  previous 
phases  of  a  guided  flight  of  an  aerodynamic  decelerator  system  that  was  developed  earlier  for  a  simple  one¬ 
dimensional  wind  model  can  be  modified  to  accommodate  more  complex  models  as  well.  Specifically,  the  paper 
presented  modifications  involving:  i)  variable  (with  altitude)  downwind  component  of  the  wind  based  on  no  prior 
knowledge  of  the  surface-layer  winds  or  incorporating  ground  wind  data  provided  by  the  target  ground  weather 
station  on-line,  ii)  variable  downwind  and  crosswind  components,  again  based  on  online  wind  estimates  at  the 
current  altitude  and/or  surface-layer  wind  magnitude  and  direction  provided  by  the  target  ground  weather  station, 
and  iii)  variable  3D  wind  uplinked  in  real  time  by  the  preceding  system.  Even  the  original  guidance  and  control 
architecture  assures  unprecedented  into-the-wind  touchdown  accuracy  of  about  10m  CEP  with  a  maximum  miss 
distance  of  30m,  demonstrated  in  over  a  hundred  drops  of  the  miniature  autonomous  parafoil  delivery  system 
Snowflake.  Modifications  presented  in  this  paper  will  allow  users  to  utilize  more  complex  tactical  scenarios,  e.g. 
intentional  landing  with  a  substantial  crosswind  component  or  operating  in  the  mountainous  areas  with  significant 
variations  in  the  vertical  component  of  the  wind,  while  preserving  the  superb  touchdown  accuracy. 
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